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Abstract 

We have developed a formalism, based on the Fourier- Bessel expansion, that facilitates 
the evaluation of matrix elements involving nucleon recoil operators, such as appear in 
serveral exotic forms of neutrinoless double beta decay (/3/3oi/)- The method is illustrated 
by applying it to the "charged" majoron model, which is one of the few that can hope 
to produce an observable effect. From our numerical computations within the QRPA 
performed for ^^Ge, ^'^Se, ^°°Mo, ^^^Te and ^^^Nd nuclei, we test the validity of approxi- 
mations made in earlier work to simplify the new matrix elements, showing that they are 
accurate to within 15%. Our new method is also suitable for computing other previously 
unevaluated PPqu nuclear matrix elements. 
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Since the classic majoron model of Gelmini and Roncadelli [|I| was excluded by the LEP 
measurement of the invisible width of the Z boson 0, ^, several new possibilities have 
been proposed for generating exotic neutrinoless double-beta decay with the emission of a 
massless majoron (/3/3m)- An interesting class of such models are those where the majoron 
is either a scalar ||^ or a vector boson that has lepton number —2 to balance that of 
the emitted electrons. These "charged majoron" models (CMM's), so-called because the 
majoron carries the unbroken f/(l) charge of lepton number, are probably the only ones 
which have a chance of producing (5(5m at a rate which could be observed in the present 
generation of experimental searches 0, |3, H- 

/^From the point of view of nuclear physics, the CMM's are interesting because their 
predicted rate of f^i^u depends on different matrix elements than in previously considered 
types of double beta decay, which had not been hitherto computed. In a previous publica- 
tion we made a first step toward computing the CMM matrix elements by neglecting a 
number of operators, an approximation which, while widely used in the literature, should 
be validated by more accurate computations. The aim of this letter is twofold. First, we 
elaborate a formalism, based on the Fourier-Bessel expansion introduced in ref. [^], that 
allows us to evaluate in a rather simple way the matrix elements of nuclear recoil opera- 
tors, such as appear in the CMM. Second, to quantify the contributions of the operators 
neglected previously, we perform the corresponding numerical calculations in the QRPA 
for ^^Ge, 82 ^e, i°0Mo, ^^STe and ^^"^Nd nuclei. 

We begin by describing the matrix elements of interest and the approximations that 
were previously made. In the CMM of ref. 0, on which we will focus in this letter, the 
inverse half-life for /3/3M-decay is 

[T(0+ ^ 0+)]-^ = ^7L^cm|-M+„ - M-^\ (1) 
where Qcm is the effective majoron coupling defined in ref. 

= /(Q - ^1 - ^^f n he.F{^M^^. (2) 



is the kinematical factor in natural units [T^ and 
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with 



OcAf (fc, rnm) = i J d^k^ ■ Yfl^^e*''-'""™, (4) 

are the nuclear form factors for charged majoron emission corresponding to the exchange 
of two neutrinos with masses M± 0. Here 



v{k; M, 



±) 



vr I Mi 
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co'±(cj± + /i) k{k + jj,) 



+ 



1 d 



(5) 



k + M| ) , and fx is the average excitation energy of the intermediate nuclear 



states. The operator 



Rnm 



Cm - CnCTm) + igvgA{(Tn X + D„ X CT^) + g1{T)r 



(6) 



with 



D„ 



1 



2M 
1 

2M 



CTn ■ (2p„ + q„), 

(2p„ + q.„ + ifwCTn X qr. 



(7) 



is the impulse non-relativistic approximation for the recoil term, obtained through the 



standard Foldy-Wouthuysen transformation [jTT|-|T^. M is the nucleon mass, fw = 4.7 
is the effective weak-magnetism coupling constant and p„ = — "iV^ and q„ = — are, 
respectively, the nuclear and lepton operators acting at the point r„, with the under- 
standing that operates only on the nucleon wave functions, while acts only on 
the exponential factor in eq. (W). Therefore one can make the replacements 



qn 



which yield^ 

iMk ■ YRnm 



-9l 



9v9a 



9l 



-{(Tn ■ k){crm ■ k) + (cr„ • k)( 

""m ■ Pm) \(^n 

■ p„)(cr„ ■ k)] 

^k ■ (cr„ X Pm + p„ X <Tm) + fw[{(^n " k) (cT^ " k) — ( 

iPn - Pm) -k+k"^ 



(9) 



^In eqs. (|^) and we have neglected the induced pseudoscalar interaction term. It is relatively small 
and can be easily embodied in eq. ^ by simply multiplying the operator ((t„ • k) ((t,„ ■ k) by the factor 
(1 + hpWo/2), where hp = gp/gA — 2M/mJ ~ 1/20, and Wq is /3/3 transition energy of the order of 
MeV. 
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In our previous work [0] we adopted the approximation made in another context by 
Doi et al, 

ik-YRnm = gA^o-n- (Tmi2gyfw + g^), (10) 

and which was also used by ref. p|. Although a series of different justifications was given 
for this procedure, and we will expand upon them below, operationally it means to neglect 
in (^ all terms containing the nucleon momenta, as well as the term g'^k'^, and to use the 
relation 

(cr„ ■ k)(cr„ ■ k) = -A;V„, ■ cr„. (11) 

However the nuclei momenta are necessarily of the same order as the momentum transfer 
k, and from this point of view the trustworthiness of the estimate ([T0|) is questionable and 
it should be checked numerically. In the remainder of this letter we shall describe a method 
for so doing, and illustrate its use by computing for several nuclei the contributions to 
the matrix elements of k ■ YRnm that were previously neglected. 
To begin, we rewrite the operator OcM{k,rnm) in the form 



OcM{k,rnm) = Ocltik, 



with 



2 

0^t(A;,r„„) = ll / d^ke'^'-^-lian ■ k)(o-„ • k) - 2(o-„ ■ k)(o-™ ■ p^)], 

OZAk. rnm) = -^-Jf j dn^e'''-'- {z2k ■ o-„ x p„ + /^[(o-„ ■ k)(o-„ ■ k) - k\a„. ■ a^)]} 

2 

Oll{k,v^J = -|^|rfn,e^k.r_^2p„-k+A;2). (13) 

In deriving the last result we used the fact that the terms containing p^ contribute in the 
same way as the analogous terms containing p„, i.e., the contribution of k ■ (cr„ x p^) is 
equal to that of — k ■ (cr^ x p„), etcl^ To obtain the corresponding matrix elements 

M^, = 2^(0;iE / k'dkv{k;M^)O-cM{k,r^m)\0t), (M) 



^Our relative signs differ in this equation and in eq. (g) from the corresponding ones of ref. because 
of a sign error in the latter. 

•^This can be easily shown from the multipole expansion that follows. 
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where |0^) and |0j ) are the initial and final nuclear states, respectively, and x stands for 
AA, Av and yy, we do the following. First, we perform the multipole expansion in ([T3|) 
and integrate over the angular coordinates dVLk. Second, we use the angular momenta 
algebra and introduce a complete set of intermediate states \Jl). In this way we get 



Mt1:^ = ^ fdkkMk;M^) 

^'^ LL'J^ 

X {25iz.'^^+^(-l)^^i-'(^l|J)(0;||TiJ(fc,r,o-)||J:)(J:||YJ(fc,r,o-■p)||0+) 
- fcLL'J-2(Ll|J)(L'l|J)(0;||Tz.j(A;,r,o-)||J:)(J:||Tiv(fc,r,o-)||0r)}, 



(15) 



I 

LL'JJi 



X {~Uk6LL'{0^\\TLAk,r,a)\\J:){J:\\TLj{k,r,a)\\0t) 

+ /^fcLL'J-2(Ll|J)(L'l|J)(0+||Tz.j(A;,r,o-)||J:)(J:||Tiv(fc,r,o-)||0n 



- 2f-l 



v^^jt i |(Ll|L')(0;||TUfc,r,^)||J:)(J:||Wfc,r,p)||Or)|, 

(16) 

-Mom" = ^ / dkkMk; M±) E 

X {-M,.,,(-l)"(0;i|Y,(A:,r,l)||J:)(J:||Y,(A:,r,l)||0+) 

- 2^^+^LJ-i(Ll|J) (0+||Yj(A;,r,l)||J:)(J:||Tz.j(fc,r,p)||0+)}. (17) 



where J = y/2J + 1, (ivl|J) is a short notation for the Clebsh-Gordon coefficient (L010| JO) 
and 

YLM{k,r,l) = jL{kr)YLM{i), 
YLM^k, r, cr ■ p) = j L{kr)YLM{r){<y ■ p), 
TLjM(fc, r,cr) = iL{kr){(T ®Yl{y))jm, 

TLjAf(/c,r,p) = jL(A;r)(p®rL(r))jM, (18) 
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are spherical one-body operators |Tl]], the reduced matrix elements of which appear in 
eqs. ( p!5| - p!7| ). The explicit dependence on the energies ojj^ of the intermediate states \Jl) 
in the latter equations can be restored, if desired, by replacing /i ojjn in eq. (^. The 
total matrix element is clearly 

and in the approximation of eq. (|10[) this simplifies to 

X {0^\\TLj{k,r,a)\\J:){J:\\TLj{k,r,CT)\\Ot). (20) 

In what follows we shall compute both the the full matrix element and this approximation 
in order to compare the two. 

The results derived so far are valid for any nuclear model; it remains only to evaluate 
the reduced matrix elements (0^| |Cij(fc, r) 1 1 J_^) and (J^| |0^j(A;, r)| |0j'') of the operators 
OLj{k,r) and (9^j(A;,r) hsted in (|T8|) . Within the QRPA formulation, after solving the 
BCS equations for the intermediate nucleus |16], the transition matrix elements become 

(Oj||CLj(fc,r)||J;) = -'^{p\\OLj{k,r)\\n) [VpUnXpn-aJ + UpVnYpn-aj] , 

pn 

{J:\\O'i^j{k,r)\\0t) = -J2{p\\0'Lj{k,r)\\n)[upVnXp^,^j + VpU^Yp^,^j], (21) 

pn 

and the BCS approximation results from the substitution: 

5:(0+||OUfc,r)||J:)(J:||OUA:,r)||0+) 

a 

^ Y,{p\\OLj{k, r) I \n) {p\\0'Lj{k, r) | \n)upVnVpUn. (22) 



pn 



For harmonic oscillator radial wave functions, the reduced single-particle pn form 



factors are JTT], |T^] 

{p\\Yj{k,r,l)\\n) = {A7v)-^W{pn;J0J)R\pn;J,k), 

{p\\TLj{k,r,CT)\\n) = {AiT)-^^W{pn; L,l, J)R''{pn; L,k), 

TLj{k,r,p)\\n) = (An)--^ \wt\pn; L, J)R^-\pn; L, k) 
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+ Wf^\pn; L, J)R^+\pn; L, k)] , 
Yj{k,r,(T ■p)\\n) = {A7T)-^wt\pn;J)R''-\pn;J,k) 
+ Wt\pn; J)R^+\pn; J, k)] , 



with the angular partsQ 

W{pn; L, S, J) 
wi^\pn;L,J) 

Wi'^^pn-J) 

and the radial parts 

R\pn;L,k) i 
R^^\pn-L,k) -- 



V2SJLininip{lnL\lp) I L S J 

^i{-i)'v+3M\jLi;j;j^{i^ + i T k)HipL\in T 1) 

1 jn In J \ \ J-n ~F 1 ^P j 



1 i i 

2 2 
jn In ~F 1 



^71 ~F 1 Jn o 



Jp 



/•OO 

R^{k;lp,np,ln,nn) = Unp,i^{r)un„,i„{r)3L{kr)r^dr 



= ± ' {^^^" + 2n„ + 2 T /p, rip, /„ T 1, 

+ (2n„ + 1 ± l)^R^{k;lp,np,ln =F l,^n ± 1)} , 



(23) 



(24) 



(25) 



where u = Muo/h is the oscillator parameter. We use here Qa = Qv, ^ value almost uni- 
versaly adopted for the effective axial vector coupling constant, both in single and double 
beta nuclear decays. As in our previous work we made the numerical calculations in 
the simplifying limit of M+ = oo, which implies Ai'^^j = 0, and M_ = 100 MeV, and we 
followed the QRPA procedure described in refs. |jT6|, |18| .0 For the discussion that follows 



*We use here the angular momentum coupHng 

^ There was a misprint in ref. [Q requiring the following replacement: 



LSJ^ 



^To assess the reliability of the nuclear model employed here, it is worthwhile to compare the recently 
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we will only need the parameters s and t, defined as the ratios between the T = 1, S = 
and T = 0, 5* = 1 coupling constants in the particle-particle (PP) channels and the pair- 
ing force constants, i.e., s = 2vPp / [vP"'"' (p) + ^^^'^(n)] and t = 2vf^ /[vP'''''{p) + <'^^^(n)]. 
(For a value of s = 1 the isospin symmetry is restored within the QRPA.) 

Table gives the contributions of the individual matrix elements in (|T2|) and (|T^ and 
their combinations for the charged majoron /3/?m decay in ^^Ge. Three different results 
are shown: i) the unperturbed or BCS values (second column), ii) the QRPA calculations 
when only the particle-hole interaction is considered, i.e., with s = t = (third column) 
and iii) the full QRPA calculations with s = 1 and t = 1.25 (fourth column). From the 
results for the Gamow- Teller like operators A;^cr„ ■ cTm/M and cr„ ■ kcr^ ■ k/M we see 
that the approximation (11) is always quite good. This is also true for the individual 
contributions, coming both from the natural parity (vr = (— )'^) and the unnatural parity 
(tt = {—)'^^^) virtual states J". The last two rows show that the approximation (10 



IS 



also quite reasonable. This is clearly because the contributions coming from the velocity 
dependent operators, as well as from the Fermi-like operator k'^/M, are relatively small. 
It should be noticed that the operators — 2(cr„ ■ k)(<Tm ■ Pm.)/M and — 2ik ■ <t„ x p^/M 
(— 2p„ ■ k/M and —k'^/M) are forbidden by the L and J selection rules for the natural 
(unnatural) parity intermediate states in eqs. ([T3|), (0) and (|T^. Moreover, while all 
intermediate states contribute coherently for the velocity independent operators, the same 
is not true of the velocity dependent operators. The latter are also more sensitive to the 
PP correlations than the former, to the extent that their contributions pass through zero 
for t < 1.25. Table H shows the comparison of the exact QRPA results from eq. (plOf) , for 



several (3(3 decaying nuclei, with those resulting from the approximation (|20|) . We have 
used here s = 1 and t = t^y^ as discussed in ref. |]l8l, i.e., t = 1.25, 1.30, 1.50, 1.40 and 
1.85 for "^^Ge, ^"^Se, ^^^Mo, ^"^^Te and ^^^Nd, respectively. In all cases the exact matrix 
elements are larger in magnitude than the approximate ones by < 15%. 

We would also like to point out that our formulation is applicable to matrix elements 



reported result Tl'j.^ = (4.3ti i[stat.] ± 1.4[syst.]) x 10^^ y for the /3/32i. half-life measurement in ''^^Ca, 
with the predicted value (using the same nuclear structure technique) T^^2 = 2.8 x 10^^ y (l8). When 
the calculation was done only the lower limit of the 2^ decay half-life, T^y2 ^ 3.6 x 10^^ y, was known 
experimentally. On the other hand, the nuclear shell model |^] restricts the corresponding half-life from 
above, T^y^ < lO^o y. 
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that appear in the usual f3f3o,^-decay, as well some supersymmetric contributions. For the 



former, the neutrino mass term has been evaluated within this framework in refs. H, 1 



In addition, the recoil matrix element ||20||-[E3 



< = (-' / (2^'' ■ ^ °- + °" ^ M*T7I) ) ■ 

which contributes to neutrinoless j3f3 decay in models where the right-handed {V + A) 
exists along with the usual left-handed {V — A) one, can be easily calculated from eq. 
(0) by making the substitution 

9.9..(*;MJ^i^^. (27) 

(The finite nucleon size effect and the short-range two-nucleon correlations can also be 
incorporated in a simple way .[]) Furthermore in supersymmetric models p4| a matrix 
element arises with a neutrino potential of the form 

hn{r) ~ - / dk-^Mkr) = f T^^^e^''^, (28) 



IT J uj{uj + jj) J 2^"^ uj{uj + jJ,) 

whose extra powers of k'^ relative to the usual neutrino potential show that this is similar 
to the recoil-like terms we evaluated above, making this matrix element also amenable to 
calculation by our method. 

In summary, we have developed a formalism especially suited for computing nuclear 
form factors that contain nuclear recoil operators, regardless of the nuclear model, and we 
applied it to the matrix elements for charged majoron emission in the QRPA. We thereby 
found that the approximations used in our previous paper to simplify the nuclear form 
factor (^) work quite well, to within a relative error of 15%. 



^ So far the matrix element M-^^'^'' has been handled as a two-body operator, which leads to rather 
complicated analytical expressions (see eqs. (2.20) to (2.27) in ref. |^ and eqs. (3.65) to (3.68) in 
ref. H). 
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Table 1: Contributions of the individual matrix elements and their combinations for the 
charged majoron decay in ''^Ge. 



Operator 


BCS 


ph-QRPA 


full-QRPA 




0.200 


0.165 


0.112 


(Tn ■ kcr^ • k/M 


0.067 


0.055 


0.037 


-k'^/M 


0.067 


0.049 


0.035 


-2(o-„-k)(o-^-p^)/M 


-0.026 


-0.014 


0.007 




0.016 


0.014 


-0.004 


-2p„ ■ k/M 


0.006 


0.006 


0.006 


~9a^ ' nCm ~ CnCr m) 


0.041 


0.041 


0.044 


-igvgA^ ■ (cr„ X Dm + D„ X cr„) 


0.640 


0.530 


0.348 


-(^^k • (D„ - Dm) 


0.073 


0.054 


0.041 


^k ■ Y/j„m 


0.754 


0.624 


0.433 


5(^Pcr„ ■ (Tm{'2gvfw + gA)/^M 


0.694 


0.571 


0.387 



Table 2: Results for the exact and approximated matrix elements, given eqs. (|19D and 
(pop, respectively. Both the particle-hole and the particle-particle channels have been 
considered in the QRPA calculations. 



Operator 


^^Ge ^^Se ^^^Mo ^^^Te ''''^Nd 


^k • ~V Rnm 

gAk'^o'n ■ (Tra{2gyU- + gA)/^M 


0.433 0.437 0.444 0.412 0.345 
0.387 0.388 0.380 0.359 0.288 
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